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A semimetal-insulator transition in the Hubbard model on the honeycomb lattice is studied by 
using the dynamical mean field theory. Electrons in the honeycomb lattice resemble the Dirac 
electron liquid and for weak interactions the system is semimetal. With increasing the local in- 
teraction a semimetal-insulator transition occurs. We find a nonanalytical structure of the phase 
transition which consists of a first-order transition line ending in a second-order transition point 
and high-temperature crossover line. A phase separation of semimetal and insulator occurs at low 
temperatures. Maxwell construction is performed to determine the first order transition line. The 
phase diagram is also presented. 
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I. INTRODUCTION 

The study of correlation driven metal-insulator transi- 
tion (MIT) has attracted great attention in recent years. 
The MIT is realized in a number of transition metal ox- 
ides and organic compounds by application of the pres- 
sure or chemical substitutions^ Far from the transition 
point the metallic phase is well described by the Fermi 
liquid theory while in the insulating phase the electrons 
are localized. When the magnetic frustration is large, 
the MIT occurs in the paramagnetic phase. The metal- 
lic state is a Fermi liquid with a renormalized mass. 
The renormalized mass increases as the transition is ap- 
proached. This is the essence of the Brinkman-Rice the- 
ory of the MIT. 2 The Brinkman-Rice scenario of the MIT 
is substantially developed by the dynamical mean field 
theory (DMFT) of the Hubbard model in the paramag- 
netic phase^i 5 - The essential features of the MIT stud- 
ied within the DMFT is a nonanalytical structure of the 
phase transition, which consists of the first-order transi- 
tion line ending in a second-order transition point and the 
high temperature crossover linei 5 -^ The nonanalytical 
structure of the MIT has been observed experimentally, 8 
as well as has been confirmed by cluster DMFT— and 
by other techniques ^ 

Recently, the experimental realization of a single layer 
of graphite, known as graphene ; 12 ' 13 has brought up re- 
newed interest in the low temperature physics of the elec- 
trons on the honeycomb lattice. In the honeycomb lattice 
at half filling the noninteracting Fermi surface collapses 
into the edge points of the Brillouin zone. The tight- 
binding dispersion exhibits the Dirac cone near these 
points, and the density of states (DOS) at the Fermi 
level vanishes. The electrons on the honeycomb lattice 
closely resemble the one of massless Dirac fermions in 
2+1 dimensions. In particular, the Hubbard model on 
the honeycomb lattice can be considered as an asymp- 
totic infrared massive quantum electron dynamics in 2 + 1 
dimensions.^ Therefore, the electrons on the honeycomb 
lattice provide a condensed matter analogy of relativistic 
physics of electrons. 



The honeycomb lattice is also a basis structure of a 
number of materials such as magnesium diboridei^ or 
layered nitride superconductor o 16 ' 17 ' 18 . Electron correla- 
tions constitute the essential properties of these materials 
such as unconventional super conductivity ! ' The emer- 
gence of electron correlations and the specific features 
of the honeycomb lattice structure underlies the mate- 
rial properties. The Hubbard model on the honeycomb 
lattice is a minimal model to describe the emergence of 
electron correlations in the specific lattice structure. 

For weak local interactions the electrons on the hon- 
eycomb lattice always stay in the paramagnetic state. 
There is no presence of superconducting or magnetic in- 
stabilities at weak coupling^ The interaction only renor- 
malizes the Fermi velocity. With increasing the local 
interaction the renormalized velocity decreases. When 
the renormalized velocity vanishes the electrons are lo- 
calized and the state is insulating i2I It is a scenario of 
the semimetal-insulator transition (SMIT) on the honey- 
comb lattice. The transition occurs in the paramagnetic 
phase, and it is a version of the Mott MIT. However, due 
to the absence of the quasiparticle mass, the Brinkman- 
Rice scenario of the MIT cannot be applied to the honey- 
comb lattice. Parallel to the SMIT, at low temperatures 
the local interaction in the honeycomb lattice can lock 
electrons with different spins into the different sublattices 
that creates a magnetic long-range order. The numeri- 
cal simulations also find a semimetal - antifcrromagnetic 
insulator transition (SMAFIT) on the honeycomb lat- 
tice at low temperaturesi 22 ' 23 ! 24 This transition is a type 
of the Slater MIT which is driven by long-range order. 
In the honeycomb lattice the Mott and the Slater MIT 
compete with each other. However, when the magnetic 
frustration is strong they can destroy the magnetic long- 
range order, and there is only the Mott transition. The 
previous DMFT studies of the SMIT on the honeycomb 
lattice showed that the SMIT is a second-order transition 
in the paramagnetic phase ; 21 ' 25 whereas the variational 
calculations showed a first-order transition characteristic 
of the SMITi^ In this paper we reexamine the SMIT on 
the honeycomb lattice by the DMFT. We find a nonan- 
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alytical structure of the SMIT on the honeycomb lattice 
which has not been pointed out in the previous DMFT 
studies . 21 ! 25 The nonanalytical structure is reminiscent to 
the phase structure of the MIT in the square or the Bethe 
lattices where the DOS at the Fermi level is finitei 5 -^ It 
also consists of the first-order transition line ending in a 
second-order transition point and the high temperature 
crossover line. However, in contrast to the square or the 
Bethe lattices, the SMIT on the honeycomb lattice does 
not accompany the appearance and the disappearance of 
a quasiparticle peak at the Fermi level. The absence of a 
quasiparticle state at the Fermi level is a specific feature 
of the electron dynamics in the honeycomb lattice. The 
SMIT on the honeycomb lattice is accompanied with the 
appearance and the disappearance of a pseudogap struc- 
ture near the Fermi level. 

The outline of the present paper is as follows. In Sec. II 
we describe the DMFT for the Hubbard model in the 
honeycomb lattice. Numerical results are presented in 
Sec. III. In Sec. IV conclusion and remarks are presented. 

II. MODEL AND METHOD 

We consider the Hubbard model on the honeycomb 
lattice. The Hamiltonian of the system is 

H = -t ^ 4a C ja~ f J '^2 c L C ia + U ^2 n ^ n il> i 1 ) 

<i,j>,a i i 

where cj CT and Ci a are the creation and the annihilation 
operator of electrons at site i with spin a. rii a — c\ a Cic, 
is the density operator, t is the hopping integral between 
nearest neighbor sites i and j. We will take t as the unit 
of energy. U is the local interaction and /i is the chemi- 
cal potential. In the following we will consider only the 
half filling case, i.e. // = U/2. The honeycomb lattice 
has a specific feature in the band structure, where the 
Fermi surface of noninteracting electrons at half filling 
is just the edge points of the Brillouin zone. The tight- 
binding dispersion near these points exhibits the Dirac 
cone like the relativistic electrons and the DOS at the 
Fermi level vanishes. The noninteracting electrons on the 
honeycomb lattice is a semimetal which has a zero gap 
at the Fermi level. The specific feature distinguishes the 
honeycomb lattice from other lattices such as the square 
or the Bethe lattices where the DOS at the Fermi level 
is finite. Apparently, in the honeycomb lattice the con- 
cept of effective mass is not appropriate to describe the 
electron properties. As a consequence, the standard de- 
scription of the Mott MIT is not valid in the honeycomb 
lattice. To reveal the nature of the SMIT in the honey- 
comb lattice we use the DMFTM^ The DMFT is exact 
in infinite dimensions. However, for two dimensional lat- 
tices the DMFT neglects nonlocal correlations. The clus- 
ter DMFT studies for a square lattice have shown that 
the key features of the MIT are already captured by the 
single-site DMFT.— Thus, one can expect a similarity 
for the honeycomb lattice. 



Since the honeycomb lattice is a Bravais lattice with a 
basis of two lattice sites, the electron Green function can 
be written in the form of matrix 2x2 

G(k,tw„) = [Go^k.iwn) - £(Mw n )]-\ (2) 

where uj n = (2n — 1)itT is the Matsubara frequency, 
S(k, iu> n ) is the self energy, and Go(k, iui n ) is the bare 
Green function of noninteracting electrons. The bare 
Green function is 

Gf'(k.^)-(^ (3) 

where e(k) = — 2t exp(ifc x /2) cos(^/Sk y /2) — texp(ik x ). 
Equation © is just the Dyson equation. Within the 
DMFT, the self energy is approximated by a local func- 
tion of frequency, i.e., 

S(k, iu n ) w 8 a pE(iw n ). (4) 

Note that within the DMFT the off diagonal elements of 
the self energy are neglected. These elements vanish in 
infinite dimensions. In finite dimension lattices they are 
indeed nonlocal correlation quantities. The self energy 
E(iu> n ) is determined from the dynamics of single-site 
electrons embedded in an effective mean field medium. 
Once the effective single-site problem is solved the self 
energy is calculated by the Dyson equation 

Y,{iu n ) = 5 _1 (iun) - G~ x {iuj n ), (5) 

where Q(u>) is the bare Green function of the effective 
single site and represents the effective mean field acting 
on the site. G(u>) is the electron Green function of the ef- 
fective single site. The self consistent condition requires 
that the Green function G(u>) of the effective single site 
must coincide with the local Green function of the origi- 
nal lattice, i.e., 

G(iuj n ) = — ^ G aa (k, ico„), (6) 

k 

where N is the number of lattice sites. Equations ©- 
([6]) form the self consistent system of equations for the 
lattice Green function and the self energy. They are prin- 
cipal equations of the DMFT. In order to solve the effec- 
tive single site problem we use the exact diagonalization 
technique! 5 -^ The exact diagonalization maps the effec- 
tive single site problem into an Anderson impurity model 

flAIM = Y^ E P b U b P<r +^2 V p( b U c ^ +h.C.) 
pa pa 

-ficlc^ + Un^ni, (7) 

where the local impurity represented by cj., c a couples to 
a bath of free conduction electrons represented by b p(J , 
b pa with dispersion E p via a hybridization V p . In the 
exact diagonalization the effective medium Green func- 
tion Q(iu) n ) is approximated by the corresponding Green 
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function G n ,(iuj n ) calculated within the Anderson impu- 
rity model of a finite number of bath levels 

™ s W I 2 

g n .(iun)- 1 = iv n + li-Y, ■ v > ( 8 ) 
^ tuj n - E p 

where n s — 1 is the number of bath levels. The model 
parameters E p , V p are determined by minimization of 
the distance function 

M 

d = M E \"n\- k \Q{iUn)- X - GnAiUnT 1 ?- (9) 
n=l 

The parameter k if chosen large (fc > 1) enhances the 
importance of the lowest Matsubara frequencies in the 
minimization procedure. In particular, we take k = 3 
in the following numerical calculations. When the model 
parameters E p , V p are obtained we solve the Anderson 
impurity model by the exact diagonalization and obtain 
the local Green function G(ioj n ), and then the self energy 
Ti(iu) n ). Thus, we obtain a closed self consistent system 
of equations for determining the electron Green function 
within the DMFT. 
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III. NUMERICAL RESULTS 

We solve the DMFT equations by iterations'^ Most 
calculations are performed with M = 1024 positive Mat- 
subara frequencies and lattice size of 64 x 64 sites. For 
very low temperature (for instance T = 0.005) we take 
M = 2048. The exact diagonalization of the Ander- 
son impurity model is performed with n s = 7. We have 
checked the agreement between Q(iu) n ) and G na {i&n) and 
found a good agreement for whole model parameter range 
under consideration. We find two typical solutions, one 
is semimetal and the other is insulator. The DOS of 
these solutions is presented in Fig. [TJ The semimetal 
solution is characterized by a pseudogap near the Fermi 
level, while the insulator solution opens a wide gap at the 
Fermi level. In the semimetallic phase the DOS shows 
two Hubbard subbands and the pseudogap structure be- 
tween them. In contrast to the square or the Bethe lat- 
tices, no Kondo quasiparticle peak appears at the Fermi 
level. In the honeycomb lattice the DOS of noninter- 
acting electrons linearly vanishes at the Fermi level, so 
that the Kondo-singlet formation resulted in the effective 
single-site problem is suppressed^ The feature of the 
noninteracting DOS is retained in the interacting case, 
so that the relativistic properties of the electrons in the 
honeycomb lattice are maintained as far as the system 
does not approach the SMIT. The local interaction only 
renormalizes the Fermi velocity which can be seen by the 
increase of the slope of the DOS near the Fermi level 
when the interaction increases. When the slope of the 
DOS near the Fermi level becomes very large, it closes 
the pseudogap and the system transforms to the insu- 
lating phase. At the point of the SMIT the pseudogap 



FIG. 1: Density of states p{lu) for various U at temperature 
T = 0.01. The first two top panels plot the semimetal solution 
while the last two panels plot the insulator solution. 

structure disappears and leaves a wide gap in the DOS. 
In the insulating phase the DOS exhibits only two Hub- 
bard subbands separated by the gap. The SMIT scenario 
in the honeycomb lattice is reminiscent to the MIT in the 
square or the Bethe lattices. The crucial different feature 
is the absence of the Kondo quasiparticle state at the 
Fermi level in the honeycomb lattice. As a consequence, 
the Brinkman-Rlce scenario of a divergence of the effec- 
tive mass at the transition point is not valid. Instead, 
in the honeycomb lattice the renormalized Fermi veloc- 
ity vanishes when the system crosses the transition point. 
This SMIT scenario was also observed in the DMFT stud- 
ies with the iterated perturbation theory as the impurity 
solver. 2 - However, at finite temperature we observe a co- 
existence of the semimetallic and the insulating solutions 
at intermediate interactions which has not been pointed 
in the previous DMFT studies i 21 i 25 

In Fig. [2] we present the imaginary part of the self 
energy for various interactions in both the semimetallic 
and the insulating phases. The slope of Im£(ia;) for oj — > 
is identical to the slope of ReY,(ui + i0 + ) for lu — > 0. The 
renormalized factor of the Fermi velocity is 



dReT,(uj + i0 + ) 
Figuref5]shows that in the semimetallic phase fcnE(iu>) 
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FIG. 2: (Color online) The imaginary part of the self energy 
for various interactions U at temperature T = 0.01. The filled 
symbols are the self energy in the semimetallic phase, while 
the open symbols are the self energy in the insulating phase. 



as lu —>■ 0. It leads the DOS to be vanished at the 
Fermi level. As the interaction increases the slope of 
Im£(iu;) for u> — > increases, so that the renormalized 
factor Z gradually decreases. When the system trans- 
forms to the insulating phase the renormalized factor Z 
vanishes. Apparently, in the insulating phase the concept 
of the Fermi velocity is not valid and the use of Eq. [TU] 
does not make sense. Nevertheless, at the SMIT point 
the slope of ImS(iw) for lu — > abruptly changes from 
a negative large value to a positive large value. It is a 
specific feature of the SMIT in the honeycomb lattice. In 
the square or the Bethe lattices, the slope of lmE(iu>) for 
io — > continuously changes as the system crosses the 
MIT£&1 

At low temperatures we have found the insulator so- 
lution for U > U c i(T), and the semimetal solution for 
U < U c2 (T). In the range U cl (T) < U < U c2 (T) both 
insulator and semimetal solutions coexist. The nonana- 
lytical structure in the honeycomb lattice is very similar 
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FIG. 3: (Color online) The quantity Q = -l/Im£(i-7rT) as a 
function of U at different temperature T. 



to the one in the square or the Bethe lattices^^ This 
suggests that the SMIT in the honeycomb lattice is a 
first order transition. In order to reveal the nonanalyti- 
cal structure of the phase transition we use the method 
proposed by Tong et al£& It is based on the observation 
that for a fixed temperature the formal dependence of a 
thermodynamical quantity Q on the interaction U is a 
multivalued function }iq(U). Usually, hQ{U) has a "Z"- 
or "S" -shaped structure. The signal of the nonanalyti- 
cal structure is the discontinuity of Q(U) in the normal 
calculations. To obtain the multivalued function /iq(J7), 
instead of Q — hq{U) we transform it to a self-consistent 
equation 

Q = h Q (U-X[A-Q}), (11) 

where A and A are parameters which are chosen so that 
Q is single valued with respect to U even if the original 
hqiJJ) is a multivalued function. Equation (fTTj) is em- 
bedded into the DMFT self-consistent equations. In the 
following calculations we take Q = — l/ImS(i7rT) which 
is the inverse of the imaginary part of the self energy 
at the first Matsubara frequency. This quantity is pro- 
portional to ImS(icJi)/wi w <9ReE(ijj)/9ijj| w _»0! which is 
a renormalized contribution to the Fermi velocity. In 
Fig. [3]we present Q as a function of U at various temper- 
atures. One can see the "Z"-shaped structure of fiQ(U) 
at low temperatures. At high temperatures }iq(U) is 
a single-value function. In the insulating phase Q van- 
ishes, while in the semimetallic phase it is finite. The 
vanishing of Q also indicates the vanishing of the renor- 
malized Fermi velocity. At low temperatures in the range 
Ud{T) < U < U C 2{T) we find an additional solution to 
the semimetal and the insulator solutions. An example 
of these three coexistent solutions are plotted in Fig. 2] 
One solution is semimetal (filled circles) and the other 
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FIG. 4: (Color online) Three coexistent solutions at U = 9 
and T = 0.01. The imaginary part of the self energy is plotted. 
The semimetal solution is presented by the blue filled circles, 
the insulator solution is presented by the red open circles, 
and the metastable solution is presented by the brown filled 
triangles. 
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FIG. 5: (Color online) The double occupation D as a function 
of U at various temperatures. The solid line, the filled cycles, 
triangles and squares are the double occupation at temper- 
ature T = 0.01, 0.05, 0.1, 0.15, respectively. The vertical 
dashed line is a Maxwell construction for T = 0.01. 

is insulator (open circles). The additional solution (filled 
triangle) is found between them. It fills the positive slope 
piece of the multivalued function }iq{U). At high fre- 
quencies the additional solution closely approaches to the 
insulator solution, while at low frequencies it behaves like 
metallic. With the three coexistent solutions the func- 
tion Iiq{U) is continuous but multivalued in the region 
U cl (T) < U < U c2 {T), as can be seen in Fig. [3 

However, not all three solutions are stable. To find a 
stable solution we compare the free energies of the three 
coexistent solutions. The free energy can be calculated 
via the double occupation^ 

F(U,T) = F(0,T)+ / D(U',T)dU', (12) 
Jo 

where F(U,T) is the free energy and D(U,T) — (ft^n^) 
is the double occupation. The double occupation can be 
calculated by the exact relation^ 

D = ~ J2 G(iiv n Micj n )e^" 0+ . (13) 

n 

The double occupation is a measure of the portion of 
lattice sites which are occupied by electrons with both 
spins and characterizes the mobility degree of electrons 
in the lattice. For U = 0, D = 0.25 for U = and 
for U — > oo, D = 0. The double occupation is often 
used to reveal the first order phase transitions^ The 
DMFT shows that within a stable phase the double oc- 
cupation decreases as U increases, and it exhibits a dis- 
continuity when the system crosses a first order transition 
line,^ Indeed, D = dF/dU and dD/dU = d 2 F/dU 2 . 
dD/dU < shows a stability of phase. The variable pair 
D — U is analogous to the inverse density and pressure 
in the conventional liquid-gas transition theory ] 28 ' 30 In 
Fig. [5] we present the double occupation D as a func- 
tion of U at various temperatures. It shows that the 



FIG. 6: Phase diagram of the SMIT in the honeycomb lattice. 



double occupation in the insulating phase is indepen- 
dent on temperature. It means that the thermal fluc- 
tuations do not affect the degree of the electron mobility 
in the insulating phase, because the gap strongly pre- 
vents the mobility of electrons. At high temperatures 
D(U) is a single-value function, while at low tempera- 
tures D(U) has "Z" -shaped structure like the quantity 
Q. The obtained function D(U) is very similar to the 
one in the square and Bethe lattices^ At high temper- 
atures the double occupation D{U) is smooth; thus, the 
system only crosses from semimetal to insulator. For 
low temperatures, at U c \ the double occupation D(U) 
is continuous, but its slope is discontinuous. It means 
that the transition at U c \ is a second-order transition. 
At U C 2 the double occupation is smooth; thus the system 
only crossovers from a semimetal to a metallic-like phase. 
However, comparing the free energy, one can see that the 
additional metallic-like phase has highest free energy, and 
therefore it is unstable. This feature is similar to the 
conventional liquid-gas transition, where the liquid and 
the gas phases coexist. The semimetallic phase is sta- 
ble for U < U C (T) 1 while the insulator one is stable for 
U > U C (T). U C (T) can be found by a Maxwell construc- 
tion Fsm{Uc,T) — Fj(U c ,T). At low temperatures as 
U crosses U C (T) from below, a stable semimetallic phase 
transforms into a stable insulating phase. This transi- 
tion is accompanied with a finite jump AD of the double 
occupation, that it is a first-order transition. The finite 
low-temperature SMIT in the honeycomb lattice is of first 
order. When T — > 0, one can expect that U c — > U C 2, and 
at U c there is no jump of the double occupation. The 
zero temperature SMIT in the honeycomb lattice is of 
second order, 2 - However, this second-order phase transi- 
tion is special. It emerges from the metastable coexistent 
phases. At zero temperature near the phase transition a 
metastable phase still exists. 

We summarize the results in the phase diagram pre- 
sented in Fig. [5J At low temperatures the semimetal- 
lic phase persists up to U C 2(T), while the insulating 
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phase exists down to U C %(T). Both lines U C \(T) and 
U C 2{T) terminate at the point Usmit- I n the region 
Ud(T) < U < U C 2(T) three phases coexist. One phase is 
scmimetallic, and the other is insulating. The third phase 
is unstable. Actually, the SMIT transition occurs along 
the line U C (T) which presents a first-order transition. 
The line U C (T) also terminates at the second-order tran- 
sition point Usmit- In the region U C \{T) < U < U C (T) 
the semimetallic phase is stable and the insulating phase 
is metastable, whereas in the region U C {T) < U < U C 2(T) 
the insulating phase is stable and the semimetallic phase 
is metastable. The phase separation is the essential fea- 
ture of the SMIT on the honeycomb lattice. This is rem- 
iniscent of the MIT in the square or the Bethe lattices. 
The feature may be considered as a common correlation 
effect regardless of the lattice structure. This may be un- 
expected because the MIT in the square or Bethe lattices 
is accompanied with the appearance and disappearance 
of a quasiparticle peak which is formed by the Kondo 
effect at the Fermi level. In the honeycomb lattice the 
Kondo effect is suspended and the Kondo quasiparticle 
peak at the Fermi level is absent. The SMIT on the 
honeycomb lattice is accompanied with appearance and 
disappearance of the pseudogap near the Fermi level. Ap- 
parently, the phase separation in the MIT is just an emer- 
gence of electron correlations in the boundary of metallic 
and insulating phases without involving a specific tran- 
sition driven mechanism. Near the zero temperature the 
first order phase transition occurs at U c /t ss 10. The pre- 
vious DMFT calculations adopting the iterated perturba- 
tion theory for single-site problem obtained U c /t = 13.3 
at T = Oi^i It is well known that the iterated perturba- 
tion theory usually overestimates the critical valued The 
DMFT calculations for infinite dimension hyperdiamond 
lattice obtained U c /t 8. 5. 25 In graphene samples^ 
U/t ~ 2 4 is far from the SMIT. One may expect that 
the graphene is well described by the Dirac liquid theory 
with a renormalized velocity. 

IV. CONCLUSION 

In this paper we study the SMIT in the honeycomb 
lattice by using the DMFT. In contrast to the square or 



Bethe lattices, the SMIT in the honeycomb lattice occurs 
without involving the appearance and the disappearance 
of a quasiparticle state at the Fermi level. It is accom- 
panied by the appearance and the disappearance of a 
pseudogap near the Fermi level. Far from the transition 
point the semimetallic phase is a Dirac electron liquid 
with a renormalized Fermi velocity, while in the insulat- 
ing phase the electrons are localized. When the system 
approaches the SMIT point from below the renormalized 
Fermi velocity vanishes. We found a nonanalytical struc- 
ture of the phase transition. It consists of a first-order 
transition line ending in a second-order transition point 
and high-temperature crossover line. At low tempera- 
tures the phase separation between the semimetallic and 
the insulating phases occurs. It suggests that the phase 
separation in a common feature of the Mott MIT regard- 
less of a specific transition driven mechanism. In two- 
dimensional lattices the DMFT neglects nonlocal correla- 
tions. The cluster DMFT calculations for a square lattice 
shows that the first order characteristic of the Mott MIT 
are already captured in the single-site DMFT»iS How- 
ever, the nonlocal correlations modify the shape of the 
transition finest One may expect the same feature in 
the honeycomb lattice. Moreover, in the honeycomb lat- 
tice the magnetic instability may compete with the Mott 
MIT at low temperatures. It requires a further study, at 
least, how strong short range magnetic correlations affect 
the Mott MIT. 
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